Ковариационная матрица регрессионных коэффициентов
2019.02.11
В курсе про регрессию пишут, что ковариационная матрица для коэффициентов регрессии нужна для рассчета стандартных ошибок оценки. Но как себе представить ковариацию коэффициентов? На самом деле, их истинных значений с нулевой вариацией мы не знаем, а знаем лишь оценки на основе данных. Т.е. сами коэффициенты являются случайными величинами, которые мы оцениваем. А посему точность их оценки может быть изучена.
library(tidyverse)
library(ggplot2)
swiss_scaled = swiss %>% mutate_all(scale)
ms = lm(Fertility ~ ., swiss_scaled)
ms %>% summary
Call:
lm(formula = Fertility ~ ., data = swiss_scaled)
Residuals:
Min 1Q Median 3Q Max
-1.22275 -0.42122 0.04029 0.32980 1.22652
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.162e-16 8.367e-02 0.000 1.00000
Agriculture -3.129e-01 1.278e-01 -2.448 0.01873 *
Examination -1.648e-01 1.621e-01 -1.016 0.31546
Education -6.704e-01 1.409e-01 -4.758 2.43e-05 ***
Catholic 3.476e-01 1.177e-01 2.953 0.00519 **
Infant.Mortality 2.511e-01 8.901e-02 2.822 0.00734 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5736 on 41 degrees of freedom
Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10vcov(ms) %>% round(3) (Intercept) Agriculture Examination Education Catholic Infant.Mortality
(Intercept) 0.007 0.000 0.000 0.000 0.000 0.000
Agriculture 0.000 0.016 0.005 0.007 -0.003 0.003
Examination 0.000 0.005 0.026 -0.013 0.011 0.000
Education 0.000 0.007 -0.013 0.020 -0.008 0.002
Catholic 0.000 -0.003 0.011 -0.008 0.014 -0.002
Infant.Mortality 0.000 0.003 0.000 0.002 -0.002 0.008# а еще посмотрим на зависимости между переменными
ms2 = lm(Fertility ~ (.)^2, swiss_scaled)
ms2 %>% summary
Call:
lm(formula = Fertility ~ (.)^2, data = swiss_scaled)
Residuals:
Min 1Q Median 3Q Max
-0.70157 -0.31115 -0.05445 0.25119 1.12881
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01212 0.18137 0.067 0.94715
Agriculture -0.30657 0.16000 -1.916 0.06462 .
Examination -0.20774 0.21217 -0.979 0.33511
Education -0.68458 0.20935 -3.270 0.00264 **
Catholic 0.16811 0.20964 0.802 0.42872
Infant.Mortality 0.20197 0.11445 1.765 0.08747 .
Agriculture:Examination 0.31001 0.19980 1.552 0.13091
Agriculture:Education 0.33321 0.26624 1.252 0.22009
Agriculture:Catholic 0.19914 0.21609 0.922 0.36387
Agriculture:Infant.Mortality 0.33732 0.15785 2.137 0.04060 *
Examination:Education 0.46164 0.22319 2.068 0.04703 *
Examination:Catholic -0.04082 0.28726 -0.142 0.88791
Examination:Infant.Mortality 0.31812 0.24009 1.325 0.19485
Education:Catholic -0.22894 0.32668 -0.701 0.48865
Education:Infant.Mortality 0.07530 0.27846 0.270 0.78863
Catholic:Infant.Mortality 0.09645 0.15724 0.613 0.54409
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5182 on 31 degrees of freedom
Multiple R-squared: 0.819, Adjusted R-squared: 0.7314
F-statistic: 9.352 on 15 and 31 DF, p-value: 1.077e-07Можно теперь нагенерить разных оценок коэффициентов регрессии из разных подвыборок методом бутстрепа и посмотреть на зависимость оценок одних коэффициентов от других. По идее, это должно идейно напоминать матрицу ковариции.
R = 1000
numbers = 1:nrow(swiss_scaled)
regressors = swiss_scaled %>% colnames()
regressors[1] = "(Intercept)"
coef_df = data.frame(R = numeric(),
Coeff = character(),
Val = numeric(), stringsAsFactors = F)
count = 1
for (i in 1:R){
m_current = lm(Fertility ~ ., swiss_scaled[sample(numbers, replace = T),]) %>% coefficients()
for (regr in regressors){
coef_df[count,"R"] = i
coef_df[count,"Coeff"] = regr
coef_df[count,"Val"] = m_current[regr]
count = count + 1
}
}
vcov(ms) %>% round(3) (Intercept) Agriculture Examination Education Catholic Infant.Mortality
(Intercept) 0.007 0.000 0.000 0.000 0.000 0.000
Agriculture 0.000 0.016 0.005 0.007 -0.003 0.003
Examination 0.000 0.005 0.026 -0.013 0.011 0.000
Education 0.000 0.007 -0.013 0.020 -0.008 0.002
Catholic 0.000 -0.003 0.011 -0.008 0.014 -0.002
Infant.Mortality 0.000 0.003 0.000 0.002 -0.002 0.008library(GGally)
coef_df_w = coef_df %>% spread(Coeff, Val) %>% select(-R)
cov(coef_df_w) %>% round(3) (Intercept) Agriculture Catholic Education Examination Infant.Mortality
(Intercept) 0.007 -0.001 0.001 0.004 -0.002 0.001
Agriculture -0.001 0.013 -0.002 0.008 0.004 0.000
Catholic 0.001 -0.002 0.011 -0.006 0.009 -0.002
Education 0.004 0.008 -0.006 0.030 -0.016 0.000
Examination -0.002 0.004 0.009 -0.016 0.026 0.001
Infant.Mortality 0.001 0.000 -0.002 0.000 0.001 0.010vcov(ms) %>% round(3) (Intercept) Agriculture Examination Education Catholic Infant.Mortality
(Intercept) 0.007 0.000 0.000 0.000 0.000 0.000
Agriculture 0.000 0.016 0.005 0.007 -0.003 0.003
Examination 0.000 0.005 0.026 -0.013 0.011 0.000
Education 0.000 0.007 -0.013 0.020 -0.008 0.002
Catholic 0.000 -0.003 0.011 -0.008 0.014 -0.002
Infant.Mortality 0.000 0.003 0.000 0.002 -0.002 0.008ggpairs(coef_df_w)